# total random sites to create
tot <- nrow(restdat)
id <- stri_rand_strings(tot, 5)
dts <- range(restdat$date)
ext <- bbox(tbpoly)
lon <- rnorm(10 * tot, mean(reststat$lon), sd(reststat$lon))
lat <- rnorm(10 * tot, mean(reststat$lat), sd(reststat$lat))
tmp <- SpatialPoints(cbind(lon, lat),
proj4string = crs(tbpoly)
) %>%
.[tbpoly, ] %>%
.[sample(1:nrow(.@coords), tot, replace = F), ]
restdat_rnd <- tibble(id) %>%
mutate(
date = rnorm(tot, mean(restdat$date), sd(restdat$date)),
date = round(date, 0),
top = sample(c('hab', 'wtr'), tot, replace = T)
)
reststat_rnd <- tibble(id) %>%
mutate(
lat = tmp$lat,
lon = tmp$lon
)
resgrp <- 'top'
restall_rnd <- left_join(restdat_rnd, reststat_rnd, by = 'id')
names(restall_rnd)[names(restall_rnd) %in% resgrp] <- 'Restoration\ngroup'
# base map
ext <- make_bbox(reststat_rnd$lon, reststat_rnd$lat, f = 0.1)
map <- get_stamenmap(ext, zoom = 10, maptype = "toner-lite")
pbase <- ggmap(map) +
theme_bw() +
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank()
)
# map by restoration type
pbase +
geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = `Restoration\ngroup`), size = 4, pch = 21)
# map by date
pbase +
geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = factor(date)), size = 4, pch = 21)
# barplot of date counts
toplo <- restall_rnd %>%
group_by(date)
ggplot(restall_rnd, aes(x = factor(date))) +
geom_bar() +
coord_flip() +
theme_bw() +
theme(
axis.title.y = element_blank()
) +
scale_y_discrete(expand = c(0, 0))
wqmtch_rnd <- get_clo(restdat_rnd, reststat_rnd, wqstat, resgrp = 'top', mtch = mtch)
save(wqmtch_rnd, file = 'data/wqmtch_rnd.RData', compress = 'xz')
head(wqmtch_rnd)
## # A tibble: 6 x 5
## stat resgrp rnk id dist
## <int> <chr> <int> <chr> <dbl>
## 1 47 hab 1 IYJz0 2960.310
## 2 47 hab 2 3OOJT 3394.566
## 3 47 hab 3 Jp5pr 3549.903
## 4 47 hab 4 GD6tr 4737.298
## 5 47 hab 5 nkm7v 4759.014
## 6 47 hab 6 Vw1c6 5440.625
##
# plots
# combine lat/lon for the plot
toplo <- wqmtch_rnd %>%
left_join(wqstat, by = 'stat') %>%
left_join(reststat_rnd, by = 'id') %>%
rename(
`Restoration\ngroup` = resgrp,
`Distance (dd)` = dist
)
# restoration project grouping column
resgrp <- 'top'
restall_rnd <- left_join(restdat_rnd, reststat_rnd, by = 'id')
names(restall_rnd)[names(restall_rnd) %in% resgrp] <- 'Restoration\ngroup'
# extent
ext <- make_bbox(wqstat$lon, wqstat$lat, f = 0.1)
map <- get_stamenmap(ext, zoom = 12, maptype = "toner-lite")
# base map
pbase <- ggmap(map) +
theme_bw() +
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank()
) +
geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = `Restoration\ngroup`), size = 4, pch = 21) +
geom_point(data = wqstat, aes(x = lon, y = lat), size = 2)
# closest
toplo1 <- filter(toplo, rnk %in% 1)
pbase +
geom_segment(data = toplo1, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)
# closest five percent
fvper <- max(toplo$rnk) %>%
`*`(0.2) %>%
ceiling
toplo2 <- filter(toplo, rnk %in% c(1:fvper))
pbase +
geom_segment(data = toplo2, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)
# closest all combo
toplo3 <- toplo
pbase +
geom_segment(data = toplo3, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)
Get weighted average of project type, treatment (before, after) of salinity for all wq station, restoration site combinations.
salchgout_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'sal', yrdf = yrdf, chgout = TRUE)
salchg_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'sal', yrdf = yrdf)
save(salchgout_rnd, file = 'data/salchgout_rnd.RData')
save(salchg_rnd, file = 'data/salchg_rnd.RData')
head(salchgout_rnd)
## # A tibble: 6 x 3
## # Groups: stat [2]
## stat cmb cval
## <int> <chr> <dbl>
## 1 6 hab_aft 25.09977
## 2 6 hab_bef 24.25426
## 3 6 wtr_aft 24.05831
## 4 6 wtr_bef 24.86349
## 5 7 hab_aft 24.82157
## 6 7 hab_bef 24.98083
head(salchg_rnd)
## # A tibble: 6 x 4
## stat hab wtr cval
## <int> <fctr> <fctr> <dbl>
## 1 6 hab_aft wtr_aft 24.57904
## 2 6 hab_aft wtr_bef 24.98163
## 3 6 hab_bef wtr_aft 24.15629
## 4 6 hab_bef wtr_bef 24.55887
## 5 7 hab_aft wtr_aft 24.62957
## 6 7 hab_aft wtr_bef 25.25801
Get conditional probability distributions for the restoration type, treatment effects, salinity as first child node in network.
wqcdt_rnd <- get_cdt(salchg_rnd, 'hab', 'wtr')
head(wqcdt_rnd)
## # A tibble: 4 x 5
## hab wtr data crv prd
## <fctr> <fctr> <list> <list> <list>
## 1 hab_aft wtr_aft <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 2 hab_aft wtr_bef <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 3 hab_bef wtr_aft <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 4 hab_bef wtr_bef <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
Discretization of salinity conditional probability distributions:
salbrk_rnd <- get_brk(wqcdt_rnd, qts = c(0.33, 0.66), 'hab', 'wtr')
salbrk_rnd
## # A tibble: 8 x 5
## hab wtr qts brk clev
## <fctr> <fctr> <dbl> <dbl> <dbl>
## 1 hab_aft wtr_aft 26.52550 0.4091514 1
## 2 hab_aft wtr_aft 30.03970 0.8643254 2
## 3 hab_aft wtr_bef 25.88986 0.3671466 1
## 4 hab_aft wtr_bef 29.78396 0.8493238 2
## 5 hab_bef wtr_aft 26.85391 0.4598276 1
## 6 hab_bef wtr_aft 30.24440 0.8849303 2
## 7 hab_bef wtr_bef 26.21855 0.4170020 1
## 8 hab_bef wtr_bef 29.98880 0.8717005 2
A plot showing the breaks:
toplo <- dplyr::select(wqcdt_rnd, -data, -crv) %>%
unnest
ggplot(toplo, aes(x = cval, y = cumest)) +
geom_line() +
geom_segment(data = salbrk_rnd, aes(x = qts, y = 0, xend = qts, yend = brk)) +
geom_segment(data = salbrk_rnd, aes(x = min(toplo$cval), y = brk, xend = qts, yend = brk)) +
facet_grid(hab ~ wtr) +
theme_bw()
Get conditional probability distributions for the restoration type, treatment effects, salinity levels, chlorophyll as second child node in network.
# get chlorophyll changes
chlchgout_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'chla', yrdf = yrdf, chgout = TRUE)
chlchg_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'chla', yrdf = yrdf)
save(chlchgout_rnd, file = 'data/chlchgout_rnd.RData')
save(chlchg_rnd, file = 'data/chlchg_rnd.RData')
# merge with salinity, bet salinity levels
salbrk_rnd <- salbrk_rnd %>%
group_by(hab, wtr) %>%
nest(.key = 'levs')
allchg_rnd <- full_join(chlchg_rnd, salchg_rnd, by = c('hab', 'wtr', 'stat')) %>%
rename(
salev = cval.y,
cval = cval.x
) %>%
group_by(hab, wtr) %>%
nest %>%
left_join(salbrk_rnd, by = c('hab', 'wtr')) %>%
mutate(
sallev = pmap(list(data, levs), function(data, levs){
# browser()
out <- data %>%
mutate(
saval = salev,
salev = cut(salev, breaks = c(-Inf, levs$qts, Inf), labels = c('lo', 'md', 'hi')),
salev = as.character(salev)
)
return(out)
})
) %>%
dplyr::select(-data, -levs) %>%
unnest
salchg_rnd <- dplyr::select(allchg_rnd, stat, hab, wtr, salev, saval)
save(salchg_rnd, file = 'data/salchg_rnd.RData', compress = 'xz')
chlcdt_rnd <- get_cdt(allchg_rnd, 'hab', 'wtr', 'salev')
save(chlcdt_rnd, file = 'data/chlcdt_rnd.RData', compress = 'xz')
chlbrk_rnd <- get_brk(chlcdt_rnd, c(0.33, 0.66), 'hab', 'wtr', 'salev')
chlbrk_rnd %>%
print(n = nrow(.))
## # A tibble: 24 x 6
## hab wtr salev qts brk clev
## <fctr> <fctr> <chr> <dbl> <dbl> <dbl>
## 1 hab_aft wtr_aft lo 10.914175 0.6595735 1
## 2 hab_aft wtr_aft lo 15.669895 0.9759047 2
## 3 hab_aft wtr_aft md 6.878474 0.5760369 1
## 4 hab_aft wtr_aft md 9.498920 0.9468640 2
## 5 hab_aft wtr_aft hi 3.584727 0.3016759 1
## 6 hab_aft wtr_aft hi 4.363590 0.7610040 2
## 7 hab_aft wtr_bef lo 9.092627 0.4355175 1
## 8 hab_aft wtr_bef lo 12.363819 0.8815996 2
## 9 hab_aft wtr_bef md 6.103452 0.3793212 1
## 10 hab_aft wtr_bef md 8.129952 0.8387278 2
## 11 hab_aft wtr_bef hi 3.473219 0.2582490 1
## 12 hab_aft wtr_bef hi 4.049840 0.7176782 2
## 13 hab_bef wtr_aft lo 11.089097 0.6925214 1
## 14 hab_bef wtr_aft lo 16.592253 0.9863017 2
## 15 hab_bef wtr_aft md 5.667113 0.3303518 1
## 16 hab_bef wtr_aft md 6.648293 0.7580916 2
## 17 hab_bef wtr_aft hi 4.077344 0.4627592 1
## 18 hab_bef wtr_aft hi 5.229716 0.8806475 2
## 19 hab_bef wtr_bef lo 9.858977 0.5375839 1
## 20 hab_bef wtr_bef lo 13.586304 0.9302304 2
## 21 hab_bef wtr_bef md 5.938513 0.4033686 1
## 22 hab_bef wtr_bef md 7.217921 0.8247996 2
## 23 hab_bef wtr_bef hi 3.616690 0.3122743 1
## 24 hab_bef wtr_bef hi 4.315727 0.7382746 2
Final combinations long format:
chlbar_rnd <- chlbrk_rnd %>%
group_by(hab, wtr, salev) %>%
nest %>%
mutate(
data = map(data, function(x){
brk <- x$brk
out <- data.frame(
lo = brk[1], md = brk[2] - brk[1], hi = 1 - brk[2]
)
return(out)
})
) %>%
unnest %>%
gather('chllev', 'chlval', lo:hi) %>%
mutate(
salev = factor(salev, levels = c('lo', 'md', 'hi')),
chllev = factor(chllev, levels = c('lo', 'md', 'hi'))
)
save(chlbar_rnd, file = 'data/chlbar_rnd.RData', compress = 'xz')
chlbar_rnd %>%
print(n = nrow(.))
## # A tibble: 36 x 5
## hab wtr salev chllev chlval
## <fctr> <fctr> <fctr> <fctr> <dbl>
## 1 hab_aft wtr_aft lo lo 0.65957353
## 2 hab_aft wtr_aft md lo 0.57603686
## 3 hab_aft wtr_aft hi lo 0.30167585
## 4 hab_aft wtr_bef lo lo 0.43551747
## 5 hab_aft wtr_bef md lo 0.37932121
## 6 hab_aft wtr_bef hi lo 0.25824904
## 7 hab_bef wtr_aft lo lo 0.69252135
## 8 hab_bef wtr_aft md lo 0.33035181
## 9 hab_bef wtr_aft hi lo 0.46275919
## 10 hab_bef wtr_bef lo lo 0.53758390
## 11 hab_bef wtr_bef md lo 0.40336858
## 12 hab_bef wtr_bef hi lo 0.31227425
## 13 hab_aft wtr_aft lo md 0.31633115
## 14 hab_aft wtr_aft md md 0.37082709
## 15 hab_aft wtr_aft hi md 0.45932817
## 16 hab_aft wtr_bef lo md 0.44608217
## 17 hab_aft wtr_bef md md 0.45940662
## 18 hab_aft wtr_bef hi md 0.45942914
## 19 hab_bef wtr_aft lo md 0.29378034
## 20 hab_bef wtr_aft md md 0.42773980
## 21 hab_bef wtr_aft hi md 0.41788835
## 22 hab_bef wtr_bef lo md 0.39264652
## 23 hab_bef wtr_bef md md 0.42143103
## 24 hab_bef wtr_bef hi md 0.42600033
## 25 hab_aft wtr_aft lo hi 0.02409532
## 26 hab_aft wtr_aft md hi 0.05313605
## 27 hab_aft wtr_aft hi hi 0.23899597
## 28 hab_aft wtr_bef lo hi 0.11840036
## 29 hab_aft wtr_bef md hi 0.16127217
## 30 hab_aft wtr_bef hi hi 0.28232182
## 31 hab_bef wtr_aft lo hi 0.01369831
## 32 hab_bef wtr_aft md hi 0.24190839
## 33 hab_bef wtr_aft hi hi 0.11935246
## 34 hab_bef wtr_bef lo hi 0.06976958
## 35 hab_bef wtr_bef md hi 0.17520039
## 36 hab_bef wtr_bef hi hi 0.26172541
Discretesize chlorophyll data, all stations:
# discretize all chl data by breaks
chlbrk_rnd <- chlbrk_rnd %>%
group_by(hab, wtr, salev) %>%
nest(.key = 'levs')
allchg_rnd <- allchg_rnd %>%
group_by(hab, wtr, salev) %>%
nest %>%
full_join(chlbrk_rnd, by = c('hab', 'wtr', 'salev')) %>%
mutate(
lev = pmap(list(data, levs), function(data, levs){
out <- data %>%
mutate(
lev = cut(cval, breaks = c(-Inf, levs$qts, Inf), labels = c('lo', 'md', 'hi')),
lev = as.character(lev)
)
return(out)
})
) %>%
dplyr::select(-data, -levs) %>%
unnest %>%
rename(
chlev = lev,
chval = cval
)
save(allchg_rnd, file = 'data/allchg_rnd.RData', compress = 'xz')
A bar plot of splits:
ggplot(chlbar_rnd, aes(x = chllev, y = chlval, group = salev, fill = salev)) +
geom_bar(stat = 'identity', position = 'dodge') +
facet_grid(hab ~ wtr) +
theme_bw()
A plot showing the breaks:
toplo <- dplyr::select(chlcdt_rnd, -data, -crv) %>%
unnest %>%
mutate(
salev = factor(salev, levels = c('lo', 'md', 'hi'))
)
chlbrk_rnd <- unnest(chlbrk_rnd)
ggplot(toplo, aes(x = cval, y = cumest, group = salev, colour = salev)) +
geom_line() +
geom_segment(data = chlbrk_rnd, aes(x = qts, y = 0, xend = qts, yend = brk)) +
geom_segment(data = chlbrk_rnd, aes(x = min(toplo$cval), y = brk, xend = qts, yend = brk)) +
facet_grid(hab ~ wtr, scales = 'free_x') +
theme_bw()